1.1 Excess mortality attributable to COVID-19

The excess mortality measure is crucial for accurately quantifying impact. The WHO defines excess mortality as “the mortality above what would be expected based on the non-crisis mortality rate in the population of interest”. Excess mortality is thus the mortality that is attributable to the crisis conditions.

To consider the impact of COVID-19, the mortality rate for year 2020, week \(t\) and age \(x\), can be represented by \[m^o_{x, t}\] and is assumed to be a result of the direct effects of COVID-19 (deaths attributable to it) and the indirect knock-on effects on health systems and society. These indirect effects impact the number of deaths due to other natural as well as external causes. The hypothetical or “counterfactual” no-COVID-19 scenario has expected mortality rate \[m^e_{x, t}\] which is estimated as the average of the most recent death rates (prior to the pandemic) or the forecast of historical rates to 2020. Excess mortality, represented by \(\mu_{x, t}\), can thus be defined as the difference: \[ \mu_{x, t} = m^o_{x, t} - m^e_{x, t} \] evaluated in rate-space to account for population changes and estimated as an absolute rate difference or percentage change. This gives a more comprehensive picture of how a location has been affected by the disease.

There is a data gap to overcome. Routine mortality data is usually received on an annual basis after the year of death or following an even longer lag. Civil registration and vital statistics (CRVS) systems differ across countries with varying timeline and quality control measures for compiling unit record cause-of-death numbers into aggregates indentified by cause, age, sex, place and period of death. In addition, differential reporting capacity and variable data quality across countries has resulted in many nations lacking the systems to provide good quality routine data in the past. Correspondingly, they also lack the capacity and data required to monitor all-cause mortality during this unprecedented pandemic. This results in a number of countries being unable to contribute to a centralized systematic mortality surveillance that would be needed to measure global, regional and country level excess mortality by the WHO.

To track excess mortality globally requires the development of tools and strategies for capturing all-cause mortality that account for the current state of irregular mortality tracking as well as imperfect country capacity. This will be achieved using three complementary exercises:

• Identifying countries that have timely and complete all-cause numbers for which the counterfactual and the observed mortality, and thus the excess, can be determined directly. This includes identifying non-conventional data sources that may provide intermediate estimates of mortality within the CRVS pipeline;
• Rapid upscale of the rate at which all-cause mortality is recorded through short term surveillance systems that are built to track all the mortality experienced during this pandemic in real time and to supplement routine systems but are not comprehensive for cause or other granular detail.
• Predicting the expected excess mortality for places with insufficient data to obtain estimates empirically. This by using statistical models that leverage the estimates from the countries where mortality is known and using them to inform plausible estimates for the countries where it is not, considering country age-, sex-, cause- and income-profiles. We note that such a solution is sub-optimal and vulnerable to the bias implicit in making inference about an unknown quantity based on incomplete information. The optimal long-term goal would be to strengthen country surveillance and to build efficient and reliable death registration systems that are at maximal coverage for all countries. In the absence of these, however, this allows for the tracking of excess mortality of more countries that can be used to generate reasonable regional and global estimates.

Over the next few sections, more detail is provided on the preliminary statistical modeling exercises.

4.1 LINEAR MODEL

The first option to predicting excess mortality is to build a log-linear model that can be applied according to a Generalized linear models (GLMs) framework. A GLM has the structure

\[g\big(\mu_i\big) = \boldsymbol{x}_i \boldsymbol{\beta}\]

where \(\mu_i\) is the expectation of response variable \(Y_i\) such that \(\mu_i \equiv \mathbb{E}\big(Y_i\big)\), \(g\) is a link function, \(\boldsymbol{x}_i\) is the \(i\)th row of the covariate matrix \(\boldsymbol{X}\), and \(\boldsymbol{\beta}\) is a vector of unknown coefficients. Akin to linear models, GLMs center around a linear predictor, \(\boldsymbol{X}\boldsymbol{\beta}\), and allow for link functions other than the identity function. GLMs also allow for any distribution from the exponential family to model the dependent variable instead of just assuming it Gaussian. We predict the log of the ratio

\[\mu_i = log(s) = log\left(\frac{Observed}{Expected + covid19}\right) = \boldsymbol{x}_i \boldsymbol{\beta}\] for \(x_i \in \{SDI, NCD, COVID\}\) where \(COVID\) represents the reported COVID mortality rate, \(NCD\) are the NCD deaths rates as described in data input as well as the interaction of the SDI and COVID. The code used to predict the scaling factor as well as the estimation of the uncertainty is provided below. Results are in a tab to follow.

#####################################################################################################################################
# Linear regression model using GLM instead of LM incase we change assumption about family
#####################################################################################################################################
rm(list = ls())
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

#####################################################################################################################################
# Load packages for all statistical analyses
pacman::p_load(highcharter, data.table, tidyr, dplyr, INLA, xgboost, mlr, countrycode)

# Can load the main data file again
load("Data/analysis_mod.Rda")

#####################################################################################################################################
# Aggregating from quarters to the sum for the entire year

df.lin      <- analysis_mod %>% 
  dplyr::select(Country, iso3, WHO_region, cov_deaths, expected, observed, Nx, ncdr, sdi) %>% 
  group_by(Country, iso3, WHO_region, Nx, ncdr, sdi) %>%
  summarize(covid = sum(cov_deaths), expected = sum(expected), 
            observed = sum(observed), .groups = "drop") %>%
  ungroup() %>% 
  mutate(excess = observed - expected, expectedc = expected + covid, 
         ratio = observed/expectedc, covidr = covid/Nx, covsdi = covidr * sdi) %>%
  arrange(iso3)

#####################################################################################################################################
# Indices to identify rows

df.lin      <- df.lin %>% mutate(row = 1:nrow(df.lin))
ind.obsl    <- df.lin %>% filter(!is.na(excess)) %>% pull(row)

#####################################################################################################################################
# Approach
# 1) Run basic linear regression model
# 2) Sample 200 sets of coefficients using distribution given by inverse of variance matrix
# 3) Generate 200 realizations of fitted scaling factors for each country
# 4) Replace scaling factor with observed for places with data
# 5) Apply scaling factors to observed covid + expected to get predicted actual
# 6) Excess is predicted actual - expected
# 7) Use samples to summarise information by country but also for regional summaries

reg.mod     <- glm(log(ratio) ~ . - 1, data = df.lin %>% dplyr::select(ratio, sdi, covidr, covsdi, ncdr))
pe          <- reg.mod$coefficients  # point estimates
vc          <- vcov(reg.mod)         # variance-covariance

set.seed(1234)
draws       <- 200

# Distribution of coefficients
simbetas    <- MASS::mvrnorm(draws, pe, vc)

# Input data to scale the coefficients
df.lin.res  <- df.lin %>% dplyr::select(dimnames(simbetas)[[2]]) %>% as.matrix()

# Distribution of fitted values (these are scaling factors from the ratio of observed to covid plus expected)
replic      <- exp(df.lin.res %*% t(simbetas))

# Use the scaling factors to obtain distribution of excess deaths
fitdist1e   <- replic*replicate(draws, df.lin$expectedc) - replicate(draws, df.lin$expected)
obsvals     <- df.lin$excess # vector of observed values

# For the non-missing, replace predictions with the actual observed
for (i in 1:draws){
  fitdist1e[ind.obsl,] <- obsvals[ind.obsl]
}

The syntax used to generate country, regional and global uncertainty bounds can be shown through the link below.

#####################################################################################################################################
# Basic functions to summarise 80% bootstrap intervals or sample from gaussian
#####################################################################################################################################

getsum <- function(x){
  tibble(mean = mean(x, na.rm = T), 
         lwr  = quantile(x, 0.1, na.rm = T), 
         uppr = quantile(x, 0.9, na.rm = T)) %>% round()
}

getsum2 <- function(x){
  c(mean(x, na.rm = T), quantile(x, 0.1, na.rm = T), quantile(x, 0.9, na.rm = T))
}

get.sample <- function(x){
  rnorm(draws, x[1], x[2])
}
#####################################################################################################################################
# Table of summarized values joined to observed and explanator data
df.lin.out   <- data.table::data.table(t(apply(fitdist1e, 1, getsum2))) %>%
  rename(y = "V1", low = "10%", high = "90%") %>%
  cbind(df.lin)  %>% 
  mutate(y    = ifelse(!is.na(excess), excess, y),
         low  = ifelse(!is.na(excess), NA, low),
         high = ifelse(!is.na(excess), NA, high), 
         source = ifelse(is.na(excess), "Predicted", "Observed")) %>% arrange(-y) %>%
  data.frame()

#####################################################################################################################################
# Indices by region

afro.indl   <- df.lin %>% filter(WHO_region == "AFRO") %>% pull(row)
emro.indl   <- df.lin %>% filter(WHO_region == "EMRO") %>% pull(row)
euro.indl   <- df.lin %>% filter(WHO_region == "EURO") %>% pull(row)
paho.indl   <- df.lin %>% filter(WHO_region == "AMRO") %>% pull(row)
searo.indl  <- df.lin %>% filter(WHO_region == "SEARO") %>% pull(row)
wpro.indl   <- df.lin %>% filter(WHO_region == "WPRO") %>% pull(row)
glob.indl   <- df.lin %>% pull(row)
#####################################################################################################################################

afro.sum   <- getsum(apply(fitdist1e[afro.indl,], 2, sum)) %>% 
  mutate(region = "AFRO", covid = sum(df.lin[afro.indl,]$covid), Nx = sum(df.lin[afro.indl,]$Nx)) 
emro.sum   <- getsum(apply(fitdist1e[emro.indl,], 2, sum)) %>% 
  mutate(region = "EMRO", covid = sum(df.lin[emro.indl,]$covid), Nx = sum(df.lin[emro.indl,]$Nx)) 
euro.sum   <- getsum(apply(fitdist1e[euro.indl,], 2, sum)) %>% 
  mutate(region = "EURO", covid = sum(df.lin[euro.indl,]$covid), Nx = sum(df.lin[euro.indl,]$Nx)) 
paho.sum   <- getsum(apply(fitdist1e[paho.indl,], 2, sum)) %>% 
  mutate(region = "PAHO", covid = sum(df.lin[paho.indl,]$covid), Nx = sum(df.lin[paho.indl,]$Nx)) 
searo.sum  <- getsum(apply(fitdist1e[searo.indl,], 2, sum)) %>% 
  mutate(region = "SEARO", covid = sum(df.lin[searo.indl,]$covid), Nx = sum(df.lin[searo.indl,]$Nx)) 
wpro.sum  <- getsum(apply(fitdist1e[wpro.indl,], 2, sum)) %>% 
  mutate(region = "WPRO", covid = sum(df.lin[wpro.indl,]$covid), Nx = sum(df.lin[wpro.indl,]$Nx))
global.sum <- getsum(apply(fitdist1e[glob.indl,], 2, sum)) %>% 
  mutate(region = "Global", covid = sum(df.lin[glob.indl,]$covid), Nx = sum(df.lin[glob.indl,]$Nx))
#####################################################################################################################################

sum.exc.df <- rbind(afro.sum, emro.sum, euro.sum, 
                    paho.sum, searo.sum, wpro.sum, global.sum) %>% 
  rename(low = lwr, high = uppr) %>%
  mutate(excess = round(mean - covid), 
         mean = round(mean), low = round(low), high = round(high),
         region = factor(region, levels = c("AFRO", "EMRO", "EURO", "PAHO", "SEARO", "WPRO", "Global")), 
         model = "Basic Linear") %>%
  dplyr::select(model, region, Nx, covid, excess, low, mean, high)

4.2 MIXED EFFECTS MODEL

A second option is to assume a linear mixed effects model to predict the excess mortality rate. In a linear mixed effects model we let \(y_j = (y_{j1}; \ldots ; y_{jn})\) denote the vector of repeated measurements for a group \(j\). We reconsider the general model \[ y_j = X_j\beta + Zjbj + \varepsilon_j\] in which \(\beta\) is a vector of population-average regression coefficients called fixed effects, and where \(b_j\) is a vector of group-specific regression coefficients. The \(b_j\) describe how the evolution of the \(j\)th group deviates from the average. The matrices \(X_j\) and \(Z_j\) are matrices of known covariates. The residual components \(\varepsilon_j\) are assumed to be independent \(N(0, \Sigma_j)\), where \(\Sigma_j\) depends on \(j\).

In this case the variable of interest is the excess mortality rate i.e., the difference between the observed and expected relative to the overall population per 100,000. As shown in the data input tab, the data are disaggregated by quarter of deaths. For fixed effects we can take the overall covid death rates for the year, the expected death rates, the sdi and quarter specific covid, NCD and communicable disease death rates as well as the quarter. Random effects assumed for the region and income group. The model, is shown in the code below with results to follow in the final tab.

#####################################################################################################################################
# Mixed effects model in INLA, similar syntax to GLM
#####################################################################################################################################

df.mixed  <- inla(excessr ~ cd_all + expectedr + sdi * covidr + ncdr + cmr + quarter - 1 + 
                  f(WHO_region, model = "iid") + f(group, model = "iid"), 
                  data = analysis_mod, 
                  control.predictor= list(compute=TRUE))

pred.in   <- cbind(df.mixed$summary.fitted.values$mean, df.mixed$summary.fitted.values$sd) 
fit.dist3 <- t(apply(pred.in, 1, get.sample)) # Assume distribution of fit is Normal with mean and sd

ind.obs <- analysis_mod %>% filter(!is.na(excess)) %>% pull(row) # indices for no data

#####################################################################################################################################
# For the observed replace predictions with actual excess rate estimates

for (i in 1:draws){
  fit.dist3[ind.obs,] <- analysis_mod$excessr[ind.obs]
}

fit.dist3e   <- analysis_mod$Nx     * fit.dist3 

#####################################################################################################################################
# Results are by quarter, need to derive draws by country and then use those to derive draws by region etc
fit.dist3ef1 <- array(dim = c(194, draws))
for (j in 1:194){
  fit.dist3ef1[j,] <- apply(fit.dist3e[(j*4 - 3):(j*4),], 2, sum)
}

#####################################################################################################################################
# Summarize with bootstrap confidence intervals based on the distribution of the draws
df.mixed.out   <- data.table::data.table(t(apply(fit.dist3ef1, 1, getsum2))) %>%
  rename(y = "V1", low = "10%", high = "90%") %>%
  cbind(df.lin) %>% filter(!is.na(y))  %>% 
  mutate(y    = ifelse(!is.na(excess), excess, y),
         low  = ifelse(!is.na(excess), NA, low),
         high = ifelse(!is.na(excess), NA, high), 
         source = ifelse(is.na(excess), "Predicted", "Observed")) %>% arrange(-y)

The code for aggregating outputs to generate region specific uncertainty can be shown using the link below.

#####################################################################################################################################
# Identify the rows for each region to aggregate the matrices accordingly

afro.ind    <- analysis_mod %>% filter(WHO_region == "AFRO") %>% pull(row)
emro.ind    <- analysis_mod %>% filter(WHO_region == "EMRO") %>% pull(row)
euro.ind    <- analysis_mod %>% filter(WHO_region == "EURO") %>% pull(row)
paho.ind    <- analysis_mod %>% filter(WHO_region == "AMRO") %>% pull(row)
searo.ind   <- analysis_mod %>% filter(WHO_region == "SEARO") %>% pull(row)
wpro.ind    <- analysis_mod %>% filter(WHO_region == "WPRO") %>% pull(row)
glob.ind    <- analysis_mod %>% pull(row)

#####################################################################################################################################
# Summarise within each aggregate, first by column (draw) and then across draws by aggregate

afro.sum.exc   <- getsum(apply(fit.dist3e[afro.ind,], 2, sum)) %>% 
  mutate(region = "AFRO", covid = sum(analysis_mod[afro.ind,]$cov_deaths), Nx = .25*sum(analysis_mod[afro.ind,]$Nx)) 
emro.sum.exc   <- getsum(apply(fit.dist3e[emro.ind,], 2, sum)) %>% 
  mutate(region = "EMRO", covid = sum(analysis_mod[emro.ind,]$cov_deaths), Nx = .25*sum(analysis_mod[emro.ind,]$Nx)) 
euro.sum.exc   <- getsum(apply(fit.dist3e[euro.ind,], 2, sum)) %>% 
  mutate(region = "EURO", covid = sum(analysis_mod[euro.ind,]$cov_deaths), Nx = .25*sum(analysis_mod[euro.ind,]$Nx)) 
paho.sum.exc   <- getsum(apply(fit.dist3e[paho.ind,], 2, sum)) %>% 
  mutate(region = "PAHO", covid = sum(analysis_mod[paho.ind,]$cov_deaths), Nx = .25*sum(analysis_mod[paho.ind,]$Nx)) 
searo.sum.exc  <- getsum(apply(fit.dist3e[searo.ind,], 2, sum)) %>% 
  mutate(region = "SEARO", covid = sum(analysis_mod[searo.ind,]$cov_deaths), Nx = .25*sum(analysis_mod[searo.ind,]$Nx)) 
wpro.sum.exc  <- getsum(apply(fit.dist3e[wpro.ind,], 2, sum)) %>% 
  mutate(region = "WPRO", covid = sum(analysis_mod[wpro.ind,]$cov_deaths), Nx = .25*sum(analysis_mod[wpro.ind,]$Nx))
global.sum.exc <- getsum(apply(fit.dist3e[glob.ind,], 2, sum)) %>% 
  mutate(region = "Global", covid = sum(analysis_mod[glob.ind,]$cov_deaths), Nx = .25*sum(analysis_mod[glob.ind,]$Nx))

#####################################################################################################################################
sum.exc.exc <- rbind(afro.sum.exc, emro.sum.exc, euro.sum.exc, 
                     paho.sum.exc, searo.sum.exc, wpro.sum.exc, global.sum.exc) %>% 
  rename(low = lwr, high = uppr) %>%
  mutate(excess = round(mean - covid), 
         mean = round(mean), low = round(low), high = round(high),
         region = factor(region, levels = c("AFRO", "EMRO", "EURO", "PAHO", "SEARO", "WPRO", "Global")), 
         model = "Mixed Effects") %>%
  dplyr::select(model, region, Nx, covid, excess, low, mean, high)

5.1 Context

For a country \(c\) with total mortality \(Y_c\) for 2020 we will report \[Y_c-E_c,\] where \(E_c\) is the expected deaths, based on a time series model applied to previous year’s data. For a country \(c\) without totaly mortality we instead use \[\widehat{Y}_c - E_c,\] where \(\widehat{Y}_c\) is the fitted value from a model based on covariates \(X_c\). In the next section we discuss models.

5.2 Models under consideration

Suppose we have the recorded number of deaths \(Y_c\) in \(2020\) and covariates \(X_c\) for each of \(C\) countries in the AMRO and EURO regions. We will explore \(3\) negative-binomial models for \(Y_c\), which take the form \[Y_c\sim \textsf{Negative-Binomial}(\mu_c, \phi),\] where \(E[Y_c]=\mu_c\) and \(\phi\) is an overdisperison parameter such that \(Var(Y_C)=\mu_c+\mu_c^2/\phi\):

• Model 1: \(\mu_c=N_c \exp(X_c\beta)\), where \(N_c\) is the population size of country \(c\),
• Model 2: \(\mu_c=E_c \exp(X_c\beta)\), where \(E_c\) is the expected number of deaths in \(2020\) for country \(c\) based on previous years’ deaths,
• Model 3: \(\mu_c=E_c + N_c \exp(X_c\beta)\).

We will then use the fitted models to predict the number of excess deaths for countries without recorded number of deaths in \(2020\).

5.3 Covariate Exploration

Before fitting the models, we explore the relationships between the covariates and the responses. First we plot \(\log(Y/N)\) vs. each of the covariates, and overlay a loess smoother:

Next we plot \(\log(Y/E)\) vs. each of the covariates, and overlay a loess smoother:

Finally we plot \(\log((Y-E)/N)\) vs. each of the covariates, and overlay a loess smoother:

5.4 Fitting the models and residual analyses

We fit each of the \(3\) candidate models in Stan using \(N(0, 5^2)\) priors on regression coefficients and \(N^+(0, 100^2)\) (i.e. a half-normal) priors on the overdispersion parameter (this prior should be explored more). Covariates were centered and scaled to have mean \(0\) and standard deviation \(1\).

We will first examine these fits through residual analyses. In particular for each model, for each of the \(C\) countries we can look at the posterior mean for \(\mu_c\) vs. the observed \(Y_c\), the posterior mean for \(Y_c-\mu_c\) vs. the posterior mean for \(\mu_c\), and the posterior mean for \((Y_c-\mu_c)/\sqrt{\mu_c+\mu_c^2/\phi}\) vs. the posterior mean for \(\mu_c\).

5.5 Leave one out metrics

We can additionally perform a leave one out cross validation residual analysis, where we fit each of the models \(C\) times, each time leaving out one of the countries, and then look at the residuals for the left out countries.

In addition to looking at residuals, using leave one out cross validation we can compute the expected log-predictive density (elpd), also known as the log CPO: \[ \sum_{c=1}^C\log\left[\int p(Y_c\mid \theta)p(\theta\mid Y_{-c})d\theta\right] \] where \(\theta\) are all model parameters, \(p(Y_c\mid \theta)\) is the predictive density for country \(c\), and \(p(\theta\mid Y_{-c})\) is the posterior conditional on all of the countries’ data except for country \(c\).

## [1] "Model 1: Offset N"

##      Estimate       SE
## elpd -685.375 15.40145
## ic   1370.750 30.80290

## [1] "Model 2: Offset E"

##       Estimate       SE
## elpd -655.3858 16.45395
## ic   1310.7716 32.90789

## [1] "Model 3: Additive"

##       Estimate       SE
## elpd -651.2184 16.01878
## ic   1302.4369 32.03756

Higher elpd is better (or equivalently lower ic, which is just -2*elpd), and so we see that models 2 and 3 have the best fits according to elpd.

5.6 Predictions

For a given country \(i\) for which we do not have the number of recorded deaths \(Y_i\) in \(2020\), we can predict \(Y_i\) using the fitted \(\mu_i\) under each of the models. For each of the countries in the AMRO and EURO regions which did not have the number of recorded deaths in \(2020\), we now plot the posterior median and \(95\%\) credible intervals for the number of excess deaths, \(\mu_i-E_i\). For the \(C\) countries for which we had the number of recorded deaths in \(2020\), we plot the observed \(Y_c-E_c\)

Here’s the same plot but slightly zoomed in:

We can also see how these excess death estimates compare to the number of reported COVID19 deaths in 2020, for the countries for which we’re making predictions.

Further we can get region level predictions of the number of excess deaths:

## [1] "Model 1: Offset N"

##   region     est   lower   upper
## 1   AMRO 1196879 1142166 1257225
## 2   EURO 1099435 1042104 1162357

## [1] "Model 1: Offset E"

##   region     est   lower   upper
## 1   AMRO 1418871 1373582 1467314
## 2   EURO 1138242 1096769 1184846

## [1] "Model 1: Additive"

##   region     est   lower   upper
## 1   AMRO 1411024 1364107 1471211
## 2   EURO 1143649 1117139 1209144

6.1 Regional comparisons across models

6.2 Regional comparisons for each model